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ABSTRACT 

Radiation hydrodynamics (RHD) simulations are used to study many astrophysical 
phenomena, however they require the use of simplified radiation transport and thermal 
prescriptions to reduce computational cost. In this paper we present a systematic study 
of the importance of microphysical processes in RHD simulations using the example 
of D-type H ll region expansion. We compare the simplest hydrogen-only models with 
those that include: ionisation of H, He, C, N, O, S and Ne, different gas metallicity, 
non-LTE metal line blanketed stellar spectral models of varying metallicity, radiation 
pressure, dust and treatment of photodissociation regions. Each of these processes 
are explicitly treated using modern numerical methods rather than parameterisation. 
In line with expectations, changes due to microphysics in either the effective number 
of ionising photons or the thermal structure of the gas lead to differences in D-type 
expansion. In general we find that more realistic calculations lead to the onset of D- 
type expansion at smaller radii and a slower subsequent expansion. Simulations of 
star forming regions using simplified microphysics are therefore likely overestimating 
the strength of radiative feedback. We find that both variations in gas metallicity 
and the inclusion of dust can affect the ionisation front evolution at the 10-20 per 
cent level over 500 kyr, which could substantially modify the results of simplified 3D 
models including feedback. Stellar metallicity, radiation pressure and the inclusion of 
photodissociation regions are all less significant effects at the 1 per cent level or less, 
rendering them of minor importance in the modelling the dynamical evolution of H II 
regions. 

Key words: stars: formation - ISM: HII regions - ISM: kinematics and dynamics - 
ISM: clouds - ISM: Bubbles - methods: numerical 


1 INTRODUCTION 


Radiative feedback has the potential to influence the mor¬ 
phological evolution of star forming regions and to induce 
or inhibit star formation (for a recent review see Dale|2015 ). 
The radiation field from massive stars ionizes the surround¬ 
ing gas, heating it and causing expansion of the hot high 
pressure region ( Spitzer||l978 Hosokawa Inutsuka|[20Q^ 
Raga et al. 2012). This can result in the dispersal of material, 

hindering the formation of stars, or can collect/destabilize 
material and potentially trigger star formation. A large num¬ 
ber of factors contribute to the effectiveness of radiative 
feedback (making models difficult, because there is a huge 
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parameter space of initial conditions and microphysical com¬ 
plexity) and feedback processes are expected to take place 
over timescales of order 1 Myr in systems of complex geom¬ 
etry (making observations difficult, due to projection effects 
and the challenge of constructing a dynamical picture from a 
single snapshot). Consequently the impact of radiative feed¬ 
back on star formation is still not clear, particularly in any 
quantitative sense (see e.g. Dale et ^|2015 ). 


Radiation hydrodynamics (RHD) calculations have 
been used to investigate the effect of radiative feedback 
in a host of “star formation” scenarios. The external ir- 


radiation of an isolated cloud (e.g. 

|Esquivel & Raga 

2007 

Gritschneder et al.|2009| Raga et al.| 

2009||Bisbas et al. 

2011 

Mackey & Lim||2011| |Tremblin et al.||20121 [Dale & Bonnel 

2012| Haworth & Harries|2012||Kinnear et al.|2014|), collect 
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and collapse (e.g. Dale et al.|2007 ) and the radiatively driven 
evolution of turbulent media (e.g. Gritschneder et aL][ 


Arthur et al.||2011 Walch et al.|2012 Tremblin et al. [2012^ 

have been studied using numerical models, providing a phe¬ 
nomenological picture of the impact of radiative feedback on 
the evolution of both the gaseous and stellar content of star 
forming regions. 

These RHD models all necessarily use (to vary¬ 
ing extents) a simplified treatment of radiation trans¬ 
port/photoionisation due to the complexity of the micro¬ 
physics and finite computational resources available. In re¬ 
cent years some effort has been invested into understanding 
how important these assumptions are. For example |Haworth| 


& Harries (2012) investigated the effect of including poly¬ 


chromatic and diffuse field radiation in models of the ra¬ 
diatively driven implosion of clouds, finding that the diffuse 
radiation field can significantly modify the results. rTremblin] 


et al. (2012) also found that the assumption of photoion¬ 


isation equilibrium can affect the results of RHD calcula¬ 
tions of the external irradiation of a turbulent medium. They 
found that non-equilibrium photoionisation was required to 
detatch the tips of elephant trunks to form Bok globules 


(Bok & Reilly 1947). More recently Sales et al. (2014) in¬ 


vestigated the effect of radiation pressure, finding that it is 
a secondary effect compared to photoionisation. |Geen et al.| 
(2015) also studied the relative impacts of winds, ionising 
radiation and supernova feedback in a series of ID models, 
which they used to summarise the energetic feedback into 
the ISM from a ISM© star. 

There are further approximations that have not been 
formally investigated, for example the assumption that the 
gas is hydrogen-only (thus neglecting cooling from forbid¬ 
den line transitions and heating and cooling from helium) 
and using a simplified thermal balance that calculates the 
temperature as a function of hydrogen ionisation fraction 
(thoug h e.g.lRaga et al.||1999| Mellema et al.|[^006| [Mackey 


Lim||2010| |Miao et ai.||2006| [Tremblin et al.||2012| com- 

pare the heating and cooling rates). There are also variations 
in the gas or stellar metallicity, dust and photodissociation 
regions/FUV heating. The impact of these approximations 
must be investigated to understand by how much and why 
simple models differ from more complex ones. 

We can easily illustrate how different approximations 
might be expected to affect the evolution of H ll regions by 
considering the classic system of a massive star at the centre 
of a uniform density medium. The star rapidly ionises a 
sphere of the surrounding gas. If the gas is hydrogen only 
and we invoke the on-the-spot (OTS) approximation, under 
which the diffuse radiation field is neglected, then the radius 
of this initial bubble is the Stromgren radius 
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where A^iy, rie and cib are the number of ionising photons 
emitted by the star, gas electron density and the case B 
recombination coefficient (for recombinations into all states 
other than the ground). The subsequent D-type expansion 
of the H II region was described by Spitzer (1978) as 


ri = rs 1 + 


7cit \ 

^rsj 
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where ci is the sound speed in the ionised gas. If we consider 


a case B recombination coefficient sensitive to temperature 
as (see equation]^ then the Stromgren radius varies 

as This will affect the expansion rate at times when 

< 1. At times when >> 1 we expect that the ex¬ 
pansion of an Hll region will vary with the temperature in 
the ionised gas as n oc and the number of ionis¬ 
ing photons as ri oc . Although these dependencies 

are quite weak, departures from simple estimates in these 
quantities may add up over time to give substantial H ll re¬ 
gion expansion differences. Furthermore, if we improve our 
model by considering gas that is not hydrogen only and in¬ 
cludes the diffuse field then both the recombination coeffi¬ 
cient and the electron density will change, further modifying 
the Stromgren radius which affects the expansion. There are 
hence multiple factors that could affect the expansion of H ii 
regions when improving the microphysical treatment. 

In this paper we aim to explore how different stellar 
and gas metallicities, radiation pressure, dust and photodis- 
socation regions affect the D-type expansion of H ll regions, 
which we interpret in terms of the simple analytic expecta¬ 
tions mentioned above. Specifically we will investigate the 
effects of these processes in galactic H ll regions rather than 
those studied in the epoch of reionisation and make the dis¬ 
tinction based on the gas densities, gas constituents and the 
ionising fluxes from the sources considered. 


2 NUMERICAL METHOD 


We use the Monte Carlo radiation transport and hydrody¬ 
namics code TORUS (Harries 2000| Kurosawa et al. 2004 


Bundle et al.|2010 Haworth et al.|2(5T2 ) to perform the cal 

culations in this paper. The TORUS RHD algorithm uses op¬ 
erator splitting to separate grid-based hydrodynamics and 
Monte Carlo photoionisation ( Haworth &: Harries|2012 ). The 
primary advantage of this method is that all of the features 
available to a dedicated Monte Carlo radiation transport 
code are available in RHD calculations. The disadvantage 
is that this approach is computationally expensive, however 
it can be efficiently parallelized using a range of techniques 


(Harries 2015). Details and testing of the RHD algorithm are 


provided in Haworth & Harries (2012) however we include 


a summary here for completeness since this paper predom¬ 
inantly explores different physical processes that the code 
can include. We also use the coupled torus-3dpdr code. 


which is discussed in detail in Bisbas et al. (2015) and also 
summarised here. 


2.1 Hydrodynamics 

We use a flux conserving, finite difference hydrodynamics al¬ 
gorithm. It is total variation diminishing (TVD) and makes 
use of the van Leer flux limiter ( van Leer|1979| and a Rhie- 
Chow interpolation scheme to prevent odd-even decoupling 
( Rhie Chow|[T983 ). TORUS is capable of treating point 
source ( Harries]|2015 ) and self-gravity, the latter of which 
is calculated using a multigrid method, though we do not 
include gravity in the models in this paper. 
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2.2 Photoionisation and thermal balance 

TORUS uses a photoionisation scheme similar to that of |Er- 


colano et ah (2003) and Wood et ah (2004) which in turn are 
based on the method presented by Lucy ( 1999| ). Packets of 
photons at constant frequency and that carry constant en¬ 
ergy e (but whose members vary in number with frequency) 
are propagated throughout the computational grid. As they 
traverse grid cells they trace a path length I and modify the 
time-averaged radiation energy density U in the cell by 


dU = 


4:71 Ju 


cAt V ^ 


( 3 ) 


where V is the cell volume, c is the speed of light, is the 
mean intensity and At is the time over which the averag¬ 
ing takes place. The update to the radiation energy density 
takes place following any photon packet ‘event’ which, as 
well as absorption, includes traversal of a cell boundary. Fol¬ 
lowing an absorption the photon packet is re-emitted with 
random frequency and direction under the principle of de¬ 
tailed balance, continuing a random walk through the grid 
until it escapes. The spectrum for diffuse field photons is 
constructed using 1000 frequencies that strategically sam¬ 
ple Lyman continuum and helium ionising photons, hydro¬ 
gen recombination lines, forbidden lines and the dust con¬ 
tinuum (the diffuse field is thus dependent on the species 
included and ionisation and temperature state of the gas). 
Once all photon packets in the calculation have escaped the 
radiation energy density is used to solve the ionisation bal¬ 
ance equation ( Osterbrock||1989 ) which, in terms of Monte 
Carlo estimators, is given by 


u(W+^) _ e 

n(X^) ~ AtVa{X^)r 


■E 


MV) 
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where n(A*), a (A*), au{X'^) and Ue are the electron num¬ 
ber density, recombination coefficient and absorption cross 
section of ion X^ and the electron density respectively. The 
approach of using the crossing of cell boundaries as a photon 
packet event has the advantage that photon packets con¬ 
tribute to the estimate of the radiation field without having 
to undergo absorption events, thus even very optically thin 
regions are properly sampled. 

For the simplest models in this paper we assume that 
the gas is either entirely atomic or ionised hydrogen. The ra¬ 
diation field is monochromatic and we use the OTS approxi¬ 
mation, including the same case B recombination coefficient 


as used by Bisbas et al. (2009) 


as = 2.7 X 10“'^ 


T 

104 


( 5 ) 


Our simple models employ a common (e.g. Gritschneder 


et al. 2009 Bisbas et al. 2011) simplified thermal bal¬ 


ance calculation where the temperature in cell j is a two- 
temperature interpolation function of the ionisation fraction 
of atomic hydrogen r/j 


Tj — Tn + rjj (Tio — T-a) 


( 6 ) 


where Tn is the prescribed temperature of fully neutral gas 
and Tio the prescribed temperature of fully ionized gas (10 
and 10^ K respectively). We retain the assumption of pho- 
toionisation equilibrium in all models. 

In contrast to the simple calculations, in the detailed 


photoionisation models we include a range of atomic con¬ 
stituents: hydrogen, helium, carbon, nitrogen, oxygen, neon 
and sulphur, for which we solve the ionisation balance using 
equation]^ The levels that we treat for metals are C (I-IV), 
N (I-IV), O (I-III), Ne (II-III), S (II-IV). The hydrogen, 
helium and CIV recombination rates used by TORUS are cal¬ 


culated based on Verner V Ferland (1996). Other radiative 


recombination rates are calculated using fits to the results 


of Nussbaumer V Storey (1983), Pequignot et al. (1991) or 
Shull V van Steenberg (1982). The photoionisation cross sec¬ 


tions of all atomic species in this paper are calculated using 


the phfit2 routine from Verner et al. (1996). We note that 


we assume rather than explicitly calculate the abundance 
of helium and metals and use this assumed abundance to 
calculate the ionisation structure. 

The detailed photoionisation model temperature is cal¬ 
culated by finding the temperature at which the heating and 
cooling rates in each cell match. The gas heating rate from 
hydrogen and helium in a given cell is calculated based on 
the sum of trajectories of photon packets through the cell 
(Wood et al. 2004). The dust heating is calculated sepa¬ 
rately, but using the same method ( Lucy||1999 ). 

The cooling rate is initially calculated for the maxi¬ 
mum and minimum allowed temperatures in the calculation 
(3 X 10^ K and 10 K respectively by default in TORUS). This is 
then refined by bisection until the cooling rate matches the 
heating rate. For gas, the cooling processes considered are 
those from free-free radiation, hydrogen and helium recom¬ 
bination and collisional excitation of hydrogen and metals. 
For dust, there is blackbody radiative cooling. We iterate 
over the ionisation and thermal balance calculations until 
the fractional change in both the ionisation fraction and 
temperature is less than 10“^. Gas and dust are thermally 
(but not dynamically) decoupled, having their thermal bal¬ 
ance solved independently. 

In the detailed calculations that are not investigating 
the effect of dust we set the dust to gas ratio to the negligibly 
low value of 10“^°. 

2.3 OSTAR2002 Spectral models from TLUSTY 

For stellar spectral models, we use the non-LTE, metal line 
blanketed, plane parallel radiation transport and hydro¬ 
static equilibrium “OSTAR2002” models of |Lanz X Huben}^ 
(2003), calculated using the code tlusty. We use three 


sets of grids for stellar metallicities of Z = 0.5,1 and 
2Z0. Each grid consists of 69 models, spanning tempera¬ 
tures from 27500 K to 55000 K and surface gravities from 
3.0 ^ log(5') ^ 4.75. The actual spectrum used in a calcula¬ 
tion is derived by interpolation between the two grid spectra 
with properties closest to that of the star in our model. In 
Figure we show an example blackbody spectum and that 
interpolated from OSTAR2002 for a star at 45000 K, a ra¬ 
dius of 10.9 R© and various stellar metallicities. The upper 
panel shows the spectrum over a large wavelength range and 
the lower panel shows the spectrum from the Lyman limit. 


2.3.1 Frequency sampling of the stellar and diffuse field 
spectra 

TORUS converts the source spectrum into a cumulative prob¬ 
ability distribution as a function of frequency. Generating a 
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Figure 1. Blackbody and tlusty spectra for a 45000 K star of 
radius 10.9 R©. The top panel shows a broad range of the spec¬ 
trum for a blackbody and star of solar metalliciity. The lower 
panel shows the ionising component of the spectrum for a black¬ 
body and stars of Z = 0.5,1 and 2 Zq. Note that the flux at the 
stellar surface is 47r times the Eddington flux. 


random number in the range 0:1 and mapping this onto the 
cumulative probability distribution then yields a frequency 
with probability appropriate to the source spectrum. The 
TLUSTY spectra consist of 193434 frequencies per model and 
the blackbody spectra consist of 1000 frequencies. The black¬ 
body spectrum is logarithmically sampled from 10 to 10^ A. 


2.5 Dust 


Historically, TORUS has been used to study discs around 
young stars where the gas and dust can be assumed to be 
thermally coupled. For modelling of H li regions, which are 
at much lower densities, we must thermally decouple the 
dust from the gas. This has been done by other codes (e.g. 


Ercolano et al.l 


2005 


Pavlyuchenkov et al. 2013), but by 


TORUS for the first time in this paper. 

We assume spherical silicate dust particles that follow 
a standard interstellar medium power law size distribution 


grains of radius a is 


Mathis et al. 1977) of the form where the number of 


( 7 ) 

where c and q are constants. 

The dust optical constants are taken from |Draine 


(1984). We use a pre-tabulated Mie-scattering phase matrix. 


At present our dust treatment does not include photoelec¬ 
tric heating or resonant line transfer. We assess the impact 
of this approximation by comparing with the more advanced 
(in terms of photoionisation) Monte Carlo photoionization 
code MOCASSIN ( [Ercolano et al.|[2003 2005) which does in¬ 
clude these dust processes. 

The model that we use for testing is the ID HII40 Lex¬ 
ington benchmark (see Eerland||l995| Ercolano et al.|[2003 


Haworth Harries||2012 ), only with the inclusion of dust. 
The dust to gas mass ratio is 10“^ and we use Draine sili¬ 
cate grains. The density is 100 mn cm“^ and the metal abun¬ 
dances used by both codes are as given in Table 

A comparison of the gas temperature distribution (all 
that really matters for these dynamic calculations, as op¬ 
posed to the line intensities) as computed by TORUS and 
MOCCASIN is given in Figure Clearly the thermally de¬ 
coupled gas and dust model from TORUS calculates a very 
similar H ii region radius and temperature to MOCASSIN. The 
extent of the ionised gas has been reduced compared to the 
simulation in which there is no dust. Beyond the ionisa¬ 
tion front there is slight heating (of order tens of Kelvin) 
in the neutral gas by TORUS, but this is weaker than in the 
MOCASSIN calculation, where there neutral gas is heated to 
~ 300 K. We do not expect this downstream heating to have 
much effect on the dynamics compared to the more dramatic 
effect on the ionisation front radius. We note that this addi¬ 
tional heating in the MOCASSIN calculation is very similar to 
photodissociation region heating, which will be investigated 
in this paper, so we will be able to gauge its impact. 


2.4 Radiation pressure 


Harries (2015) discusses and tests the treatment of radiation 


pressure by TORUS in detail. In summary the radiation pres¬ 
sure force is calculated using Monte Carlo estimators. This 
estimate is calculated as photon packets are propagated over 
the grid in the photoionisation component of the calculation 
and so if already doing the photoionisation, is essentially ob¬ 
tained for no additional computational cost. This technique 
accounts for polychromatic radiation and anisotropic scat¬ 
tering and works in both the free-streaming and optically 
thick regimes. The calculated radiation pressure appears as 
an additional force term in the hydrodynamic component of 
the calculation. 


2.6 TORUS-3DPDR 


We recently coupled TORUS with the 3D-PDR code of jBisba^ 
et al. (2012), the details and testing of which are given thor¬ 


oughly in Bisbas et al. (2015). This approach is novel since 


the UV field used in the PDR calculation is determined ex¬ 
plicitly by the Monte Carlo radiation transport (typically 
the Draine field is a free parameter for PDR codes). The 
coupled code is fully integrated, with the PDR calculations 
taking place on the TORUS grid. This means that we can 
calculate the ionisation state of multiple species in the H ll 
region as well as the chemical and thermal structure of the 
PDR and then use the resulting thermal properties as pres- 
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Figure 2. The temperature distribution in a uniform medium 
about an ionising source calculated by TORUS both with (green 
dashed line) and without dust (blue dotted line). Included is also 
the result including dust computed by MOCASSIN (red crosses). 
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Table 1. Hll region expansion model parameters. 


Variable (Unit) Value Description 


Teff(K) 45000 

nfiRQ) 10.9 



63.8 

logio(3) 

4.17 

p (rriH cm“^) 

100 

P {iTiH cm“®) 

500 

logio(He/H) 

-1 

logio(C/H) 

-3.66 

logio(N/H) 

-4.40 

logio(0/H) 

-3.48 

logio(Ne/H) 

-4.30 

logio(S/H) 

-5.05 

d/g 

1 X 10-2 

^min 

0.005 

^max 

0.25 

q 

3.3 

L (cm) 

4.4x10^9 

^cells 

256 


Source effective temperature 
Source radius 
Source mass 
Source surface gravity 
Low density model density 
High density model density 
Base helium abundance 
Base carbon abundance {zg = 1) 
Base nitrogen abundance (zg = 1) 
Base oxygen abundance (zg = 1) 
Base neon abundance (zg = 1) 
Base sulphur abundance {zg = 1) 
Dust to gas mass ratio 
Minimum dust grain size 
Maximum dust grain size 
Dust power law index 
Computational domain size 
Number of grid cells 


sure terms in the hydrodynamics. Motoyama et al. (2015) 
recently presented a 2D hydrochemical code capable of mod¬ 
elling PDR chemistry and evolving the gas dynamics. The 
difference with our code is that we can directly calculate 
the UV field over all space and can also calculate the prop¬ 
erties in the H ll region. Unfortunately such a calculation is 
computationally very expensive. 

TORUS-3DPDR uses the most recent UMIST 2012 chem¬ 
ical network database, consisting of 215 species and more 
than 3000 reactions (McElroy et al. 2013). However at 


present we use a reduced network of 33 species and 330 re¬ 
actions, primarily to reduce computational cost. Neverthe¬ 
less, RHD simulations including even this reduced network 
of species and reactions, in conjunction with photoionisa¬ 
tion, are novel. Note that we assume equilibrium in both 
the H II region and the PDR. 

Modelling of the PDR should give rise to higher temper¬ 


atures (a few hundred Kelvin rather than 10) at the bound¬ 
ary of the Hll region. This will allow us to investigate the 
effect of heating seen by dust just beyond the ionisation front 
by MOCASSIN that is not replicated by TORUS (see Figure]^. 


3 TEST MODEL SPECIFICATIONS 


The model that we use as a testbed is the classic expansion 
of an H II region about an ionising star in a uniform density 
medium (Spitzer 1978 Dyson <V Williams||1980 Hosokawa 


Inutsuka|2006 Raga et al.|2012 Bisbas et al.|2015). This 

system is ID spherically symmetric, reducing computational 
cost and thus allowing us to incorporate a large range of 
microphysics in this paper. Specifically we only consider the 
early phase expansion ( Bisbas et al.||2015 ). 

We consider a star of effective temperature 45000 K, 
radius 10.9 R© and mass 63.8 M©. These stellar parame¬ 


ters are taken from Diaz-Miller et al. (1998) and is in the 


regime where non-LTE model atmospheres are required. For 
a model star with these parameters we consider the expan¬ 
sion of the H II region into two different ambient densities of 
100 and 500 mn cm“^ which we refer to as low and high den¬ 
sity respectively. For this scenario we run a host of models 
with different physical prescriptions: simple photoionisation 
and full photoionisation (see section 2.2), non-LTE spectral 
models, different gas metallicities, radiation pressure, dust 
and treatment of photodissociation regions. We study the 
gas and stellar metallicity effects separately so that we can 
determine which, if either, the dynamics is most sensitive 
to. In reality the gas and stellar metallicities will not be in¬ 
dependent. A summary of the physical parameters in the 
model is given in Table 


4 RESULTS AND DISCUSSION 

We compare the models by tracking the ionisation front po¬ 
sition (defined as the point at which the hydrogen ionisation 
fraction is a half) as a function of time, this is shown for all 
low and high density models in Figures and We now dis¬ 
cuss each set of approximations/physical processes in turn. 


4.1 Simple versus full photoionisation 

The simple photoionisation scheme assumes an ionised gas 
temperature of 10^ K. If the temperature resulting from the 
detailed photoionisation differs from this then the expansion 
rate of the H ll region will differ, since the Stromgren radius 
varies as oc (for a case B recombination coefficient) 

and at later times during the expansion varies as n oc 
(see equation]^. The simple scheme also considers hydrogen 
only gas, which means that the electron density compared to 
a calculation with helium and metals may also differ. Finally 
the simplified scheme also uses the OTS approximation and 
so there may be differences between the case-B recombina¬ 
tion coefficient used (see equation and the net effective 
recombination coefficient for the calculation with multiple 
species. 

The top left panels of Figures and show the ionisa¬ 
tion front position as a function of time for the Spitzer an¬ 
alytic solution (equation]^ as well as for our simplified and 
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Figure 3. The ionisation front position as a function of time for all of the 100 mn cm“^ models in this paper. The top left panel compares 
the most simplified model with standard Monte Carlo photoionisation plus hydrodynamics. The top right panel includes detailed spectral 
models. Middle left varies the gas metallicity, middle right includes radiation pressure, bottom left includes dust and the bottom right 
includes PDR treatment. 
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Figure 4. The ionisation front position as a function of time for all of the 500 mn cm“^ models in this paper. The top left panel compares 
the most simplified model with standard Monte Carlo photoionisation plus hydrodynamics. The top right panel includes detailed spectral 
models. Middle left varies the gas metallicity, middle right includes radiation pressure, bottom left includes dust and the bottom right 
includes PDR treatment. 
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full photoionisation calculations with TORUS for the 100 and 
500 niH cm“^ density calculations. We also include the ana¬ 
lytic expression describing the expansion given by |Hosokawa| 


& Inutsuka (2006) (labelled H-I), which is derived by solv¬ 


ing the equation of motion of the shell of material swept up 
by the expanding H ii region 


d{MR) 

dt 


= {Pi - Po) 


( 8 ) 


which, under the assumption that Po ^ Pi, results in 


ri = Ts 



7v^c/t 

4V3r, 


4/7 


(9) 


The difference between the H-I and Spitzer solutions is that 
the former solves the equation of motion of the shell (thin 
shell approximation) and the latter considers pressure bal¬ 
ance between the Hll region and ambient medium at each 
point in time). 


4-PI Initial Stromgren radii 

Looking at the top left panels of Figures and the sim¬ 
plified models agree with the Stromgren solution for the 
position of the onset of D-type expansion in both density 
regimes, however for the detailed model this radius is smaller 
by 7.2 and 6.5 per cent in the 100 and 500 mn cm“^ models 
respectively. Given that the number of ionising photons is 
the same in each model, the difference can only come from a 
difference in the electron density or recombination rate from 
equation ^ 

We compared the electron density for the two models 
at the onset of D-type expansion and found they are very 
consistent. The difference must therefore be arising from a 
higher effective recombination rate across all species rela¬ 
tive to the case B recombination coefficient used to calcu¬ 
late the Stromgren radius analytically. The recombination 
rate (equation]^ is temperature dependent, so the difference 
could be explained if the simple and detailed photoionisa¬ 
tion schemes result in different gas temperatures. We plot 
the gas temperature as a function of radius from the ionis¬ 
ing source in Figure The simplified calculation results in 
higher average gas temperatures, implying a lower recom¬ 
bination rate and hence a larger Stromgren radius (as we 
observe in our models). 

In Figure which we will discuss more below, we 
compare our simplified and detailed photoionisation results 
against analytic solutions using gas temperatures represen¬ 
tative of the model average in the ionised gas. When we do 
this the initial I-front position of the detailed photoionisa¬ 
tion models is in good agreement with the Stromgren equa¬ 
tion. So we can conclude that the difference in Stromgren 
radius is explained by a difference in the gas temperature 
and hence recombination rate. 

We note that the Hii region radius calculated by TORUS 
for detailed photoionisation models has been shown to agree 
with other photoionisation codes, e.g. [Haworth Harris 
(2012| where we compare with the CLOUDY code (Ferland 


et al. 1998). 


Table 2. Ionizing fluxes for the different stellar spectral models. 


Stellar Model 

Metallicity 

ionising photons (xlO'^^/s) 

Blackbody 

- 

2.69 

TLUSTY 

0.5 

2.60 

TLUSTY 

1 

2.63 

TLUSTY 

2 

2.70 


4-1-2 Expansion rates 

We remind the reader that the top left panels of Figures 
and^show the ionisation front position as a function of time 
of the models being discussed presently. The StarBenclQ 
code comparison project D-type expansion test (which uses 
the simplified photoionization physics and towards which 
TORUS contributed results) found that numerical codes agree 
with the Spitzer solution early on, but eventually overtake it 
to agree more closely with the Hosokawa-Inutsuka solution. 
This behaviour is reproduced in the simple photoionization 
models in this paper, though the expansion is slightly slower 
than expected in the lower density regime. We expect that 
this is due to the comparatively low resolution (256 cells) 
used for the simulations in this paper compared to the 1024 
cells used in the ID simulations of Bisbas et al. ( 2015| ). Al¬ 
though we could increase the resolution for the simplified 
calculations (indeed our contributions to the StarBench D- 
type test were at higher resolution) this would make consis¬ 
tent resolution calculations including the PDR very compu¬ 
tationally expensive. 

The effect of moving to full photoionisation is not only 


to reduce the radius of D-type onset (see section 4.1.1), but 
also to reduce the expansion rate of the Hll region. In Fig¬ 
ure we show the temperature of the detailed and simple 
models just before the onset of D-type expansion and after 
250 kyr. Throughout most of the Hll region the tempera¬ 
ture in the ionised gas of the full photoionisation model is 
around 2000 K lower than that of the simplified model, ex¬ 
cept for the temperature peak near to the ionisation front. 
These cooler average temperatures are due to the efficient 
forbidden line cooling, which is only included in the detailed 
model. The peak in temperature close to the ionisation front 
arises where coolants with higher ionisation potentials re¬ 
combine. The simple and detailed models have average gas 
temperature of 10^ and 8000 K respectively. In Figure we 
use these average gas temperatures in the Spitzer solution 
and compare with our numerical results. This demonstrates 
that differences in the D-type expansion between the simple 
and detailed photoionisation models are simply explained by 
the differences in the ionised gas temperature. 


4.2 Stellar spectral models 

Although the different model spectra shown in Figure [^ap¬ 
pear significantly different, the integrated ionising flux over 
the spectrum is actually very similar in each case for stars 
of this effective temperature (see table in this paper and 
Figure 15 from Lanz fc Hubeny|2003| . The expansion of H ii 
regions about stars with these different spectra is therefore 


^ https://www.astro.uni-bonn.de/sb-ii/ 
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Figure 5. The radial temperature structure of the simple aud 
detailed photoiouisatiou models at the simulatiou start time aud 
after 250 kyr. The upper aud lower pauels are for the 100 aud 
500mHcm“^ ambieut deusity models respectively. 


very similar, as we show in the top right panels of Figures 
and ^ The difference in the initial Stromgren radius is neg¬ 
ligible and after 500 kyr of expansion the difference in Hii 
region radius is only of order 1-3 per cent in both density 
regimes. 

The lower panel of Figure shows that harder pho¬ 
tons are produced in the tlusty models, however this has 
very little effect on the ionisation state or dynamical evo¬ 
lution, only slightly broadening the ionisation front. This 
is in agreement with what was found by [Haworth Har-j 
(2012) in 3D models of radiatively driven implosion. 


where the effect of including harder radiation was isolated. 
Although hard radiation can have important consequences 
at galactic or cosmological scales (e.g. Shapiro et al.|[2004 


Iliev et aL]|2006 ), on scales of up to tens of parsecs we find 


that its effect is negligible. 


4.3 Gas metallicity 

The gas metallicity determines the amount of metal line 
cooling and therefore the temperature in the ionised gas. 
Since the Stromgren radius varies as T°'^^ and at later times 
ri oc the expansion rate of H ll regions is expected to 


Figure 6. Demonstrating that the differences between the sim¬ 
plified and detailed photoiouisatiou models can be explained by 
their different ionised gas temperatures. The upper and lower 
panels show results for the 100 and 500 mn cm“^ density models 
respectively. The analytic solutions are just the Spitzer equation 
with an ionised gas temperature of 10000 K and 8000 K for the 
simple and detailed photoionisation models respectively. 


be metallicity dependent. We find that this is the case in 
our models, which we show in the middle left panels of Fig¬ 
ures [ 3 ] and [ 4 ] Included are calculations at the base metal¬ 
licity (which we call zg = 1, that is the same as that used 
in the photoionisation HII40 Lexington benchmark) as well 
as zg — 0.5 (representative of earlier universe, LMC) and 
zg = 2 (representative of the Galactic centre). The lower 
metallicity expansion is much faster than the higher metal¬ 
licity, giving a zg — 0.5, 500 kyr Hll region radius larger 
than that for zg = 2 by 18 and 16 per cent in the 100 and 
500 niH cm“^ models respectively. 

In Figurej^we show the radial temperature profile of the 
different metallicity H ll regions at a time close to the onset 
of D-type expansion. The temperature in the ionised gas 
and the Stromgren radius are both decreasing as a function 
of gas metallicity. This is because at lower metallicity there 
is weaker metal line cooling. Figure shows the ionisation 
front position as a function of time for the models of different 
gas metallicity, as well as analytic plots which are generated 
from the Spitzer solution and Stromgren radius, only setting 
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Figure 7. The gas temperatures in the lOOmHcm”^ model at 
the onset of D-type expansion for gasses of different metallicity: 
0.5 (red), 1 (green) and 2 (blue) times the HII40 Lexington bench¬ 
mark metallicity. 


the ionised gas temperature to the average seen in Figure 
The agreement is good enough for us to conclude that the 
differences in H ll region expansion rate with metallicity are 
dominated by differences in the ionised gas temperature. 

Since there is good agreement between the theoretical 
expressions and simulations when the correct temperature is 
used, we fit our temperature as a function of metallicity to 
obtain an approximate, simplified thermal calculation that 
accounts for the gas metallicity, improving on the simple 
prescription given by equation 


T = Tn+ [1.1 X 10"^ - 3.8 X \&{zlzo - 0.5)° ®®® - T„] r\ 

( 10 ) 

where zjzo is the gas metallicity relative to that of the 
Lexington benchmark (a standard Milky Way star forming 
region abundance), is the neutral gas temperature and 
rj is the hydrogen ionisation fraction. An example potential 
application of this approximate simplified temperature 
calculation would be to use it 3D models such as those of 


Dale et al. (2013) or Walch et al. (2013) (both of which 
use equation to study how radiative feedback affects 
molecular clouds in regions of different metallicity, or as a 
function of cosmic time. 


Q. 

c 
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c 
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c 
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Figure 8. The numerical simulations compared with analytical 
results with fitted gas temperatures. The upper and lower panels 
are for the 100 and 500mHcm“^ models respectively. 


secondary to photoionisation in D-type expansion models. 
Further note that Raga (2015) finds that hundreds of O stars 
are required for radiation pressure to become important to 


the system dynamics, such as in the Tarantula nebula in the 
LMC. 


4.5 Dust 


In reality the gas and stellar metallicities will not be in¬ 
dependent. We investigate them separately in this paper to 
identify which of the two components dominates any differ¬ 
ences in the H ll region evolution. Based on our simulations 
we find that feedback dynamics are much more sensitive to 
the gas metallicty than changes in the stellar metallicity. 


4.4 Radiation pressure 


The ionisation front evolution for models that include radi¬ 
ation pressure is given in the middle right panels of Figures 
|3]and|4] The radiation pressure simulations are just the “de¬ 
tailed” models with the addition of the radiation pressure 
formulation developed by |Harri^ (2015). In both density 
regimes the effect of including radiation pressure is essen¬ 


tially zero (<1 per cent). This is in agreement with Sales 
et al. ( 2014| , who find that radiation pressure effects are 


Dust absorption reduces the ionising photon budget and 
hence the size of the Stromgren radius (see equation 
In Figure we also saw that it increases the temperature 
downstream of the ionisation front in the MOCASSIN results 
(though this is missing from our treatment of the dust in this 
paper, the consequences of this downstream heating effect 
will be gauged in our discussion of PDR heating in section 


4.6). The ionisation front expansion for the dust models is 
given in the bottom left panels of Figures and As we 
found when testing in Figureinclusion of dust leads to the 
initial Stromgren radius being reduced. This affects the ex¬ 
pansion rate of the H ll region at times when the Stromgren 
term is important, i.e. when < 1. Although the ionising 
photon budget is reduced, the gas temperature remains very 
similar (see Figure]^ and so at later times when n oc 
the expansion behaves in a similar manner either with or 
without dust. The effect of dust, minus any effect from down¬ 
stream heating, is therefore simply to reduce the radius at 
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which D-type expansion begins, resulting in a smaller Hii 
compared to models without dust. 


4.6 Photodissociation regions 


Models including treatment of PDRs are the most compu¬ 
tationally intensive in this paper. With these models we can 
follow ionised hydrogen, through to the atomic and molec¬ 
ular transitions (and can do the same for, e.g. carbon and 
CO and many other species). The PDR does not affect the 
temperature in the ionised gas, but leads to heating of up 
to around a few hundred Kelvin in the neutral gas out to 
relatively large distances (potentially many parsecs). 

In the derivation of Spitzer and H-I equations, there 
is a term involving the difference in the pressures in the 
ionised and neutral gas Pi — Po, where Po is dropped assum¬ 
ing Po « Pi. If the PDR heating raises the temperature 
enough then this assumption will no longer be appropriate 


and the expansion will be slower. In Bisbas et al. (2015) the 


full form of the Spitzer and Hosokawa-Inutsuka expressions 
are presented without neglecting the external pressure. 


1 dr I _ / _ /TjTo f 

a dt \ri J fioTi \ri 


- 3/4 


( 11 ) 


and 


1 dr I 
Ci dt 


\ 


Ar, 


3/2 


3r 


3/2 


l^oTi 


( 12 ) 


respectively. These expressions need to be evaluated numer¬ 
ically. In each case it is the jiiTo/jJioTi term (which is gener¬ 
ally expected to be small) that is dropped to arrive at the 
Spitzer or Hosokawa-Inutsuka results. PDR (or dust) heat¬ 
ing may change this though. 

The evolution of the H ll region radius for these PDR- 
RHD models is given in the bottom right panels of Figures 
|^and[^ In both density regimes the effect of the PDR is to 
marginally slow the H11 region expansion. The PDR heated 
gas close to the ionisation front is ~ 300 K in both density 
regimes. Substituting this temperature and an ionised gas 
temperature of 8200 K into equation pT] and comparing with 
equation!^ gives an expected difference at 500 kyr of about 1 
per cent in the Hll region extent. This is in agreement with 
the difference found in our simulations. Given the small dif¬ 
ference from downstream heating in the PDR model, we also 
expect the similar heating from dust (see Figure]^ to have 
a negligible effect. 


4.7 Consequences for RHD modelling 

We have considered the effect of many different microphysi¬ 
cal processes that are not normally included in RHD simula¬ 
tions because they make the calculation much more expen¬ 
sive (for the models including PDR treatment, by an order 
of magnitude or more). Studying all of these processes at 
once in a systematic way has the added value that they will 
not be studied across multiple papers and we can immedi¬ 
ately see which processes are important and in what ways 
current simplihed models are deficient. 


4 . 7.1 The state of simplified models 

Compared to the simplified models in this paper, every other 
process (except lowering the gas metallicity) results in a 
smaller H ll region. It is therefore possible that RHD simula¬ 
tions of more complex geometrical systems (such as turbu¬ 
lent star forming regions) using simplified microphysics are 
overestimating the power of the radiation field. In the con¬ 
text of, for example, radiative feedback and triggered star 
formation, this could mean that the ability of radiative feed¬ 
back to compress clouds and trigger star formation might 
be reduced. Conversely its ability to disperse clouds, halting 
star formation is also reduced. 

The metallicity is also particularly important, with ap¬ 
proximately a 20 per cent difference in H ll region radius af¬ 
ter 500kyr between models with metallicities 0.5 and 2 times 
our base metallicity. Radiative feedback may therefore play 
a substantially different role in the earlier universe or LMC 
compared to somewhere like the Galactic centre. Equation 
|10| a simple thermal prescription which accounts for the gas 
metallicity, could be used to investigate this in 3D models. 


4 . 7.2 Hierarehy of proeesses when modelling D-type 
expansion of galaetie Hll regions 

We have shown here that if only the dynamical evolution of 
the system needs to be modelled (as opposed to chemistry or 
synthetic observations) there is not necessarily much value in 
spending time implementing, for example, detailed spectral 
models or PDR treatment, when the effect on the dynami¬ 
cal evolution of the system is negligible. We have therefore 
constructed a simple hierarchy of processes that we find to 
be important for modelling radiative feedback in star form¬ 
ing regions, which we show in FigureOf course the effects 
will not be negligible in every scenario (for example radiation 
pressure is important early on in the formation of massive 
stars). 

The hierarchy consists of 4 tiers. Tier 1 is the standard 
basic RHD (photoionisation + hydrodynamics) code. Tier 2 
includes processes that alter the H ll region extent at 500 kyr 
by around 10 per cent or greater. Tier 2 includes processes 
that alter the Hii region extent at 500 kyr by of around 1 
per cent and tier 4 processes are untested. If improving the 
microphysics in a photoionisation + hydrodynamics code for 
non-cosmological applications this should offer a useful ref¬ 
erence as to which processes are most important to include. 
We note that in reality the gas and stellar metallicities are 
linked. Our hierarchy does not suggest that they are decou¬ 
pled, rather that if the overall metallicity varies it is the 
changes in the gas properties rather than the stellar spec¬ 
trum that affect changes in the dynamical evolution of the 
system. 


4.8 Untested features 

There are some processes that we identify that might also 
affect the D-type expansion of galactic Hll regions that we 
do not treat in this paper. 
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Figure 9. A hierarchy of importance of physical processes for 
the D-type expansion model eonsidered in this paper. There will 
be other scenarios where some of these components become sig¬ 
nificantly more important, e.g. massive star clusters such as the 
tarantula nebula in the LMC where radiation pressure plays an 
important role. Note that the percentage difference quoted is only 
after 500 kyr of H ll region evolution. 


J1..8.I X-rays and heavier metals 


Treatment of X-ray photons and a wider range of metal 
species might affect the ionisation and temperature struc¬ 
ture in the H ii region. Such treatment has been implemented 
in the Monte Carlo photoionisation code MOCCASIN (Er- 
colano et al.||20Q8 ) but is not currently available in TORUS. 
Inclusion of X-rays may lead to increased temperatures in 
the Hll region since atoms/ions with higher ionisation po¬ 
tentials will be ionised. This would lead to an increase in 
the expansion rate as we found with the lower gas metallic¬ 
ity models. 


4-8.2 Non-equilibrium photoionisation 


Including non-equilibrium photoionisaiton is not expected 
to modify the results of the calculations in this paper. It 
would allow us to capture the R-type expansion of the 
ionised gas that occurs prior to the D-type expansion that 
we study here, however R-type expansion occurs rapidly 
and with little disruption to the density field. Although 
we do not expect this assumption to affect the resulting 


density field the results of Tremblin et al. (2012) show that 


it will be important to test the effects of non-equilibrium 
photoionisation on shadowed regions. 


5 SUMMARY AND CONCLUSIONS 

We have used the Monte Carlo radiation transport and 
hydrodynamics code TORUS and its PDR counterpart 


torus-Sdpdr to study the effects of different microphysics 
on the D-type expansion of Hli regions. We have com¬ 
pared analytic solutions with the simplest RHD models, 
calculations including multiple species, the diffuse field, 
polychromatic radiation, detailed spectral models, different 
gas metallicity, radiation pressure, dust and PDR’s. Each of 
these processes has been treated directly rather than using 
a simple parameterisation. We draw the following main 
conclusions from this work: 

1) The critical factor that affects the Hll region dynamics 

when different microphysics is treated is the gas tempera¬ 
ture, which affects the Stromgren radius as Vs oc and 

the later time expansion rate as rj oc The Stromgren 

radius directly affects the expansion rate at times up until 

7 ^ . 

4rs 

2 ) Compared to simplified (hydrogen only, monochromatic, 
on the spot approximation) models, cooler average gas 
temperatures resulting from full Monte Carlo photoionisa¬ 
tion can change the extent of the Hli region after 500 kyr 
by about 10 per cent (see conclusion 1). This implies 
that simplified RHD models are overestimating the power 
of radiative feedback, be it to induce star formation or 
disperse gas in star forming regions. 

3) Varying the gas metallicity changes the gas amount 
of metal line cooling and hence the gas temperature. 
Changing the metallicity by a factor of 4 can also lead to 
a 10-20 per cent difference in the Hll region extent after 
500 kyr. The Spitzer D-type expansion solution matches 
our metallicity dependent models well if it uses the average 
gas temperatures from our simulation. We fit our models to 
produce a simplified thermal prescription that can be used 
in 3D models to incorporate heating effects of multi-species 
gas at varying metallicity. This is at negligible additional 
computational cost. It will also allow us to test how deficient 
simplified schemes are in more complex 3D applications. 

3) Dust absorption reduces the ionising photon budget, 
meaning that the onset of D-type expansion occurs at 
smaller radii (since Vs oc The Hll region temperature 

is similar, so the later time expansion rate (when the 
Stromgren radius term becomes less important) is similar 
with or without dust. 


4) Using detailed stellar spectral models does not affect 
the simulation much (of order 1 per cent difference in 
the ionised gas extent at 500kyr), since the ionising flux 
remains approximately constant and the effect of any 
harder radiation is just to smear out the ionisation front. 
Radiation pressure is also found to be dominated by the 
effects of photoionisation, in agreement with previous 
studies such as Sales et al. (2014). PDR heating of the 
downstream medium only marginally affects the H ll region 
expansion (by ~ 1 per cent at 500 kyr), in agreement with 
the difference expected from the analytic equations. 


6 ) We develop a hierarchy of processes for modelling the 
D-type expansion of H ll regions, where additional physics 
beyond the most simple model is graded based upon its 
effect upon the expansion. This offers guidance for where 
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to focus in future development of numerical models. 
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